##########################################################################################

library(dplyr)
library(data.table)
library(optparse)
library(data.table)
library(ggplot2)
library(ggrepel)

##########################################################################################

option_list <- list(
    make_option(c("--score1_file"), type = "character") ,
    make_option(c("--score2_file"), type = "character") ,
    make_option(c("--gene1_file"), type = "character") ,
    make_option(c("--gene2_file"), type = "character") ,
    make_option(c("--cgen_file"), type = "character") ,
    make_option(c("--info_multi_file"), type = "character") ,
    make_option(c("--info_TCGA_file"), type = "character") ,
    make_option(c("--type"), type = "character") ,
    make_option(c("--cnv_type"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    score1_file <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/copynumber/SCNA_unpaired/merge/CIN_IGC_Deletion.Driver_Gscore.csv"
    score2_file <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/copynumber/SCNA_unpaired/merge/CIN_DGC_Deletion.Driver_Gscore.csv"
    gene1_file <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/copynumber/gistic2/CIN_IGC_merge/all_data_by_genes.txt"
    gene2_file <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/copynumber/gistic2/CIN_DGC_merge/all_data_by_genes.txt"
    cgen_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/SCNA_driver/SCNA_driver_STAD_Deletion_candidate_Gene.csv"
    type <- "merge"
    cnv_type <- "Deletion"
    info_TCGA_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.Race.tsv"
    info_multi_file <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/copynumber/Combine_CIN_sampleInfo.tsv"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

score1_file <- opt$score1_file
score2_file <- opt$score2_file
gene1_file <- opt$gene1_file
gene2_file <- opt$gene2_file
cgen_file <- opt$cgen_file
info_multi_file <- opt$info_multi_file
info_TCGA_file <- opt$info_TCGA_file
type <- opt$type
cnv_type <- opt$cnv_type
out_path <- opt$out_path

dir.create(out_path , recursive = T)
setwd(out_path)

###########################################################################################

pri <- read.csv(score2_file)
met <- read.csv(score1_file)
cgen <- read.csv(cgen_file , header = F)

dat_gene1 <- data.frame(fread(gene1_file , header = T))
dat_gene2 <- data.frame(fread(gene2_file , header = T))
colnames(dat_gene1)[1] <- "Gene_Symbol"
colnames(dat_gene2)[1] <- "Gene_Symbol"

dat_info_multi <- data.frame(fread(info_multi_file))
dat_info_TCGA <- data.frame(fread(info_TCGA_file))

dat_info_TCGA <- subset(dat_info_TCGA,dat_info_TCGA$From=="TCGA" & dat_info_TCGA$Molecular.subtype=="CIN")
dat_info_TCGA$use_sample <- dat_info_TCGA$Tumor
dat_info_TCGA$use_sample <- gsub("-",".",dat_info_TCGA$use_sample)
dat_info_TCGA$Normal <- dat_info_TCGA$Tumor

dat_info_multi$use_sample <- paste0(dat_info_multi$Tumor,"_",dat_info_multi$Normal)
dat_info_multi$From <- "multi"

dat_info <- rbind(dat_info_multi[,c("use_sample","Normal","From")],dat_info_TCGA[,c("use_sample","Normal","From")])


diff_t <- 0.3
q_t <- 0.01

###########################################################################################
## 提取对于lauren分型的样本
if(type == "merge"){
    dat_info <- dat_info
}else{
    dat_info <- subset(dat_info , From == type)
}

use_sample <- dat_info$use_sample
dat_gene1 <- dat_gene1[,colnames(dat_gene1) %in% c( colnames(dat_gene1)[1:3] , use_sample)]
dat_gene2 <- dat_gene2[,colnames(dat_gene2) %in% c( colnames(dat_gene2)[1:3] , use_sample)]

#dat_info$use_sample <- gsub( "[-]" , "." , paste0(dat_info$Tumor , "_" , dat_info$Normal))

###########################################################################################

colnames(met)[4:5] <- c('log10Q.met','G.score.met')
colnames(pri)[4:5] <- c('log10Q.pri','G.score.pri')

if(cnv_type == "Amplification"){
    dat <- merge(met[,c('ID','symbol','chr.amp','start.amp','end.amp','log10Q.met','G.score.met')], pri[,c('ID','log10Q.pri','G.score.pri')], by='ID')
    dat$pos <- paste0( dat$chr.amp , ":" , dat$start.amp , "-" , dat$end.amp )
}else if(cnv_type == "Deletion"){
    dat <- merge(met[,c('ID','symbol','chr.del','start.del','end.del','log10Q.met','G.score.met')], pri[,c('ID','log10Q.pri','G.score.pri')], by='ID')
    dat$pos <- paste0( dat$chr.del , ":" , dat$start.del , "-" , dat$end.del )
}

## 只关注已报道
dat$symbol <- gsub( " " , "" , dat$symbol )
dat <- subset(dat, symbol %in% cgen$V1)

## 注释基因位置
dat <- merge( dat , dat_gene1[,c("Gene_Symbol" , "Cytoband")] , by.x = "symbol" , by.y = "Gene_Symbol" )
dat$show_text <- paste( dat$symbol , "\n" , dat$Cytoband )
dat$show_text <- gsub( " " , "" , dat$show_text)


###########################################################################################

dat$Gscore_diff <- dat$G.score.met - dat$G.score.pri

sig <- subset(dat, log10Q.met > -log10(q_t) & Gscore_diff > diff_t )
dat$note <- ifelse(dat$log10Q.met > -log10(q_t) & dat$Gscore_diff > diff_t, 'sig', 'non.sig')
dat$log10Q.met <- ifelse(dat$log10Q.met==Inf | dat$log10Q.met>50, 50, dat$log10Q.met)

## 显著的基因
show_symbol <- sig$symbol

#sig$symbol <- paste("italic('",sig$symbol,"')" , " " , sig$pos ,sep='')
#sig$show_text <- paste("italic('",sig$show_text,"')" , sep='')
x_lim_t <- max(c(abs(max(dat$Gscore_diff)),abs(min(dat$Gscore_diff)))) + 0.5

p <- ggplot(data = dat, aes(x = Gscore_diff, y = log10Q.met, color=note))+
     geom_point(size=0.5) +
     theme_classic()  +
     scale_color_manual(values = rev(c("#CC0000", "#000000"))) +
     xlab('GISTIC score difference') + 
     ylab(expression(paste(-log[10]," (FDR)"))) + 
     geom_hline(yintercept=-log10(q_t),color='grey50',linetype="dashed") +
     geom_vline(xintercept=diff_t,color='grey50',linetype="dashed") +
     theme(legend.position="none",
     axis.text.x = element_text(color='black'),
     axis.text.y = element_text(color='black'),
     axis.ticks = element_line(color='black')) +
     xlim(-x_lim_t,x_lim_t) + ylim(0,50) +
     # 添加标签：
     geom_text_repel(data = sig,
         #min.segment.length = Inf,
        max.overlaps = getOption("ggrepel.max.overlaps", default = 100),
        aes(label = show_text),
        color = '#CC0000',
        point.padding = 0.3,
        face = "italic" ,
        #parse = TRUE,
        size = 3) +
     #theme(text = element_text(face = "italic")) +
    labs(title = paste0(cnv_type," (",type,") "),size=2)

out_file <- paste0( out_path , "/" , "CIN_" , cnv_type,"_",type , ".Driver.pdf" )
ggsave(out_file,p,width=4.5, height=3.5)

dat$show_text <- paste( dat$symbol , "_" , dat$Cytoband )
out_file <- paste0( out_path , "/CIN_" , cnv_type,"_",type , ".Driver.tsv" )
write.table( dat , out_file , row.names = F , sep = "\t" , quote = F )

###########################################################################################
## 柱状图比较拷贝数
trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

col <- c(
    rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
    rgb(red=2,green=100,blue=190,alpha=255,max=255) 
    )

names(col) <- c("IGC" , "DGC" )

###########################################################################################

if(length(show_symbol) > 0){
    dat_gene1 <- subset( dat_gene1 , Gene_Symbol %in% show_symbol )
    dat_gene2 <- subset( dat_gene2 , Gene_Symbol %in% show_symbol )

    rownames(dat_gene1) <- dat_gene1$Gene_Symbol
    rownames(dat_gene2) <- dat_gene2$Gene_Symbol

    result <- c()
    for(gene in show_symbol){
        print(gene)
        tmp1 <- dat_gene1[gene,]
        tmp2 <- dat_gene2[gene,]

        #tmp1 <- data.frame( copy_number = 2^(as.numeric(tmp1[,-c(1:3)]) * 2) , type = "Positive" )
        #tmp2 <- data.frame( copy_number = 2^(as.numeric(tmp2[,-c(1:3)]) * 2) , type = "Negative" )

        tmp1 <- data.frame( use_sample = colnames(tmp1[,-c(1:3)]) , copy_number = as.numeric(tmp1[,-c(1:3)]) , type = "IGC" )
        tmp2 <- data.frame( use_sample = colnames(tmp2[,-c(1:3)]) , copy_number = as.numeric(tmp2[,-c(1:3)]) , type = "DGC" )

        ## 同一患者存在多个的，取中位数
        tmp1 <- merge( tmp1 , dat_info[,c("Normal" , "use_sample")] , by = "use_sample" )
        tmp2 <- merge( tmp2 , dat_info[,c("Normal" , "use_sample")] , by = "use_sample" )
        tmp1 <- tmp1 %>%
        group_by( Normal , type  ) %>%
        summarize( copy_number = median(copy_number) )
        tmp2 <- tmp2 %>%
        group_by( Normal , type  ) %>%
        summarize( copy_number = median(copy_number) )    

        p <- wilcox.test( tmp1$copy_number , tmp2$copy_number )$p.value
        #p <- t.test( tmp1$copy_number , tmp2$copy_number )$p.value

        if( p < 0.01 ){
            p_text <- trans(p)
        }else{
            p_text <- paste0( "P == " , round(as.numeric(p) , 2) ) 
        }

        tmp_res <- rbind( tmp1 , tmp2 )
        tmp_res$Gene_Symbol <- gene
        tmp_res$p <- ""
        tmp_res$p <- p
        tmp_res$p_text <- ""
        tmp_res$p_text[1] <- p_text
        result <- rbind( result , tmp_res )
    }

    ###########################################################################################
    ## 样本数量
    #dat_plot_tmp_use$id <- paste0( dat_plot_tmp_use$From , "_" , dat_plot_tmp_use$Class )
    sample_lauren_num <- result %>%
    group_by(type) %>%
    summarize( nums_lauren = length(unique(Normal)) )

    result <- merge( result , sample_lauren_num , by = "type")
    result$lauren_use <- paste0( result$type , "\n" , "(" , result$nums_lauren , ")" )
    result$lauren_use <- factor( result$lauren_use , levels = unique(result$lauren_use)[order(unique(result$lauren_use) , decreasing=T)] , order = T )

    y_max <- max(result$copy_number) * 1.5
    y_lab <- "log2(copy number/2)"
    result$type <- factor( result$type , levels = c("IGC" , "DGC") )
    ## 基于的顺序按照p值由小到达
    result$Gene_Symbol <- factor( result$Gene_Symbol , levels = unique(result[order(result$p),"Gene_Symbol"]))

    ## 注释基因位置
    result <- merge( result , dat_gene1[,c("Gene_Symbol" , "Cytoband")] , by = "Gene_Symbol" )
    result$show_text <- paste( result$Gene_Symbol , "\n" , result$Cytoband )
    result$show_text <- gsub( " " , "" , result$show_text)
    result$show_text <- factor( result$show_text , levels = unique(result[order(result$p),"show_text"]))

    sample_lauren_num <- result %>%
    group_by(type) %>%
    summarize( nums_lauren = length(unique(Normal)) )

    plot <- ggplot(result, aes(x = lauren_use , y = copy_number, color = type , fill = type)) +
            #geom_boxplot(size = 1.2 , outlier.shape = NA ) + ## 去除散点，加粗线
           # geom_jitter(position = position_jitterdodge(0.8) , size = 1) + 
            geom_violin(trim=FALSE) +
            geom_boxplot(width=0.2,position=position_dodge(0.9),fill="white",color="black")+ #绘制箱线图
            #scale_y_log10() +
            facet_grid( .~ show_text , scales = "free_x" ) +
            scale_color_manual(values=col) +
            scale_fill_manual(values=col) +
            geom_text(aes(label=p_text , y = y_max , x = 1.5),parse = TRUE,size=5 , color = "black", face='bold') +
            xlab(NULL) +
            ylab(y_lab)+
            theme_bw() +
            theme(
                legend.position = 'none',
                legend.title = element_blank() ,
                panel.grid.major=element_blank(),
                panel.grid.minor=element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 12,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 12,color="black",face='bold') ,
                axis.ticks.length = unit(0.2, "cm") ,
                strip.text.x = element_text(size = 15, colour = "black",face='italic') ,
                axis.line = element_line(size = 0.5)) 

    out_name <- paste0( out_path , "/" , "CIN_" , cnv_type,"_",type , ".Driver.copy_number.pdf" )
    if(length(show_symbol)==1){
        ggsave(file=out_name,plot=plot,width=3,height=5.4/1.5)
    }else{
        ggsave(file=out_name,plot=plot,width=length(show_symbol) * 2,height=5.4/1.5)
    }

    outcompare <- result %>%
    group_by( Gene_Symbol , type,p ) %>%
    summarize( copy_number = median(copy_number) )

    out_name <- paste0( out_path , "/" , "CIN_" , cnv_type,"_",type , ".Driver.copy_number.tsv" )
    write.table( outcompare , out_name , row.names = F , quote = F , sep = "\t" )
}